clear

% 1. Simulation settings
% Vector of seeds that creates datasets
R           =  1000;     % Number of MC simulations
seedvec     = (1:100000/R:100000);

% 1.2. TRUE PARAMETER VALUES
% Size of the datasets
nmkt        = 50;                               % number of markets = (# of cities)*(# of quarters)
prods       = 25;                               % number of brands per market
nseg        = 2;                                % number of segments in each market
nobs        = nmkt*prods;                       % number of observations
meanc       = [0 0];                            % mean of the generated vectors:
% the first is the segment latent variable, the second is the product attribute

%% Settings

% True model parameter values
betatrue    = [-1 ; -2 ; -3];               % true mean tastes on constant, dummy for segment, X1
rc_true     = [0 ; 1];                      % true random coefficient
SigsegTrue  = 0.3;

%% Set correlation: parameter that determines whether the segment is informative about the continuous characteristic or not
corrDstarX1  = 0.9;

%% This changes the numerosity of the segments
cutoff       = 1;

stdevDstar   = 1;                                % std. dev. of dstar (does not matter)
stdevX1      = 1;                                % std. dev. of X1

Chelp1       = [stdevDstar^2, corrDstarX1*stdevX1*stdevDstar; corrDstarX1*stdevX1*stdevDstar, stdevX1^2];    %varcovar matrix

rng default
ksi     = randn(nobs,1);
ksi     = ksi-mean(ksi);                        % rescale to have E[ksi]=0
const   = ones(nobs,1);

NonLinTheta = [rc_true;SigsegTrue];

% Number of parameters
nlin        = size(betatrue,1);                 % number of linear parameters
nrc         = size(rc_true,1);                  % number of random coeficients - Random coefficient Logit
nrcNL       = size(NonLinTheta,1);              % number of random coeficients - Random Coefficient Nested Logit

% Individual taste,integration of market shares using Quadrature Rule
% Note: when using quadrature we have no variation across markets
% In other words there is no randomness when integrating marketshares we
% use fixed nodes and weights as given by Heiss and Winschel (2007)

% Use now two dimensions in the simulation because I have to estimate in RC
% II two random coefficients, one on the continuous and the other one on
% the discrete characteristic
% Using two dimensions instead of one should not change anything in terms
% of precision
[nodedraws, qweight]    = nwspgr('KPN', nrc, 7);
Nnodes                  = length(qweight);
qweight                 = qweight';
qweightrprods           = reshape((qweight(ones(prods,1),:)),prods,1,1,Nnodes);

% 2. Start Simulation

% Store results
% The order of the coefficients is:
%[Constant;Segment Coefficient; Mean value of X1 (Product Attribute); Nesting parameter; Random Coefficient]
% 5=number of coefficients; 1=number of columns; 4=number of estimators (Log,NL,RC,RCNL);R=number of draws

Coefftrue            = [betatrue;SigsegTrue;rc_true];                % true parameters of the model
Coefficients         = zeros(nlin+nrcNL,1,5,R);
StandardErrors       = zeros(nlin+nrcNL,1,5,R);
ComputationTime      = zeros(1,1,4,R);

OwnElast             = zeros(nseg,nmkt,5,R);
CrossElastSameSeg    = zeros(nseg,nmkt,5,R);
CrossElastDiffSeg    = zeros(nseg,nmkt,5,R);
OneOwnElast          = zeros(1,nmkt,5,R);
OneCrossElastSameSeg = zeros(1,nmkt,5,R);
OneCrossElastDiffSeg = zeros(1,nmkt,5,R);

DiversionSameSeg    = zeros(nseg,nmkt,5,R);
DiversionDiffSeg    = zeros(nseg,nmkt,5,R);
OneDiversionSameSeg = zeros(1,nmkt,5,R);
OneDiversionDiffSeg = zeros(1,nmkt,5,R);

BiasOwnElast         = zeros(nseg,nmkt,4,R);
BiasCrossElastSameSeg= zeros(nseg,nmkt,4,R);
BiasCrossElastDiffSeg= zeros(nseg,nmkt,4,R);

FunctionValue        = zeros(1,1,5,R);
ExitFlag             = zeros(1,1,5,R);
AIC                  = zeros(1,1,5,R);
BIC                  = zeros(1,1,5,R);
% For the outside good I save the minimum, the average and the maximum for each simulation
OutsideShare         = zeros(3,R);
SaveSegment          = zeros(nobs,R);
SaveX1               = zeros(nobs,R);

for MC=1:R
    rng(seedvec(MC))
    
    % 3. Simulate data
    disp(['Simulation Numb: ',num2str(MC)])
    % Select the Seed
    
    % Generate X1 and Dstar according to a multivariate normal
    X1Dstar      = mvnrnd(meanc,Chelp1,nobs);
    
    Dstar        = X1Dstar(:,1);
    X1           = exp(X1Dstar(:,2));
    % The lognormal produces very extreme values sometimes. This
    % truncates but the distribution is not changed
    X1           = X1.*(X1 < 4) + (randn(nobs,1) + 3) .* (X1 > 4);
    multiX1      = reshape(X1,[prods,1,nmkt]);
    
    segment      = Dstar>cutoff;
    multisegment = reshape(segment,[prods,1,nmkt]);
    dumseg       = [Dstar<cutoff,Dstar>cutoff];
    multidumseg  = reshape(dumseg',[nseg,prods,nmkt]);
    multidumseg  = permute(multidumseg,[2 1 3]);
    multidumseg  = repmat(multidumseg,[1 1 1 Nnodes]);
    
    SaveSegment(:,MC) = segment;
    SaveX1(:,MC)      = X1;
    
    % random coefficient on SEGMENT and X1
    xvuRC         = X1*(nodedraws(:,1))';
    xvuSEG        = segment*(nodedraws(:,2))';

    multixvuLSEG  = reshape(xvuSEG,[prods,1,nmkt,Nnodes]);
    multixvuSEG   = bsxfun(@times,multixvuLSEG,multidumseg);
    
    multixvuLRC   = reshape(xvuRC,[prods,1,nmkt,Nnodes]);
    multixvuRC    = bsxfun(@times,multixvuLRC,multidumseg);
    
    muTrue        = xvuSEG*rc_true(1) + xvuRC*rc_true(2);
    
    Xexo          = [const segment X1];
    deltaTrue     = Xexo*betatrue+ksi;
    
    % 4. Data Structure

    % Calculate the True/Observed Market shares
    Data.nmkt           = nmkt;
    Data.prods          = prods;
    Data.nobs           = nobs;
    Data.SigsegTrue     = SigsegTrue;
    Data.nseg           = nseg;
    Data.Nnodes         = Nnodes;
    Data.muTrue         = muTrue;
    Data.dumseg         = dumseg;
    Data.qweightrprods  = qweightrprods;
    Data.xvuRC          = xvuRC;
    Data.xvuSEG         = xvuSEG;
    Data.multidumseg    = multidumseg;
    Data.multixvuRC     = multixvuRC;
    Data.multixvuLRC    = multixvuLRC;
    Data.multixvuSEG    = multixvuSEG;
    Data.multixvuLSEG   = multixvuLSEG;
    Data.nlin           = nlin;
    Data.nrc            = nrc;
    Data.nrcNL          = nrcNL;
    Data.multiX1        = multiX1;
    Data.multisegment   = multisegment;
    
    % True Market Shares
    [~, ~, sj, sjg, s0, ~, ~, ~] = TrueShare(deltaTrue,Data);
    
    % 2D reshape
    sj                         = reshape(sum(sj,2),[prods*nmkt 1]);
    sjg                        = reshape(sum(sjg,2),[prods*nmkt 1]);
    s0                         = reshape(s0,[prods*nmkt 1]);
    y                          = log(sj) - log(s0);             % log-odds ratio
    lnsjg                      = log(sjg);
    lnsj                       = log(sj);
    
    OutsideShare(:,MC)         = [min(s0);mean(s0);max(s0)]; % save outside good shares (min, mean, max)
    
    Data.logobsshare           = lnsj;
    
    % 5 Optimal instruments
    
    % Compute optimal instruments
    % Delta without ksi
    PseudoDelta                 = Xexo*betatrue;
    chamberlin                  = jacobNL(NonLinTheta,PseudoDelta,Data);
    Z                           = [Xexo chamberlin(:,2) chamberlin(:,3)]; % only retain non-zero instruments

    % 6 Estimation
    
    %% Logit
    % Weighting Matrix
    W           = (Z'*Z)\eye(size(Z,2));
    BLog        = ((Xexo'*Z*W*Z'*Xexo)\eye(size(Xexo,2)))*(Xexo'*Z*W*Z'*y);
    % Standard errors
    est         = y-Xexo*BLog;
    dgf         = (size(Xexo,1)-size(Xexo,2));
    ser         = (est'*est)./dgf;
    sst         = inv(Xexo'*Xexo);
    sterrLogit  = sqrt(ser*diag(sst));
    
    deltL       = y;
    alpha       = BLog(nlin);
    theta       = 0;
    % Logit elasticity wrt X1
    [ElaLogit, DiversionLogit]    = ElastLogit(alpha,theta,deltL,multiX1,Data);
    
    [MeanOwnElastLog,OneOwnElastLog,MeanCrossElastSameSegLog,OneCrossElastSameSegLog,MeanCrossElastDiffSegLog,OneCrossElastDiffSegLog,...
     MeanDivLog,MeanDivDiffLog,OneDivLog,OneDivDiffLog] = sumElast(ElaLogit,DiversionLogit,Data);
    
    % Function value
    FctValLog   = (est' * Z) * W * (Z' * est);
    % (size(Z),2) = number of moment conditions
    % size(BLog,1)= number of parameters
    % AIC = ln(nobs)*FunctValue-2*(numb moment conditions - numb parameters)
    AICLogit    = nobs * FctValLog - 2 * (size(Z,2)-size(BLog,1));
    % Add the two coefficients that are not estimated in the Logit (Nesting
    % parameter; Random Coefficient)
    % BIC = ln(nobs)*FunctValue-(numb moment conditions - numb parameters)*log(number of observations)
    BICLogit    = nobs * FctValLog - (size(Z,2)-size(BLog,1)) * log(nobs);
    
    % Store the results
    Coefficients(:,:,1,MC)          = [BLog;nan;nan;nan];
    StandardErrors(:,:,1,MC)        = [sterrLogit;nan;nan;nan];
    ComputationTime(:,:,1,MC)       = nan;
    OwnElast(:,:,1,MC)              = MeanOwnElastLog;
    
    CrossElastSameSeg(:,:,1,MC)     = MeanCrossElastSameSegLog;
 
    CrossElastDiffSeg(:,:,1,MC)     = MeanCrossElastDiffSegLog;

    OneOwnElast(:,:,1,MC)           = OneOwnElastLog;
    OneCrossElastSameSeg(:,:,1,MC)  = OneCrossElastSameSegLog;
    OneCrossElastDiffSeg(:,:,1,MC)  = OneCrossElastDiffSegLog;
       
    DiversionSameSeg(:,:,1,MC)     =  MeanDivLog;
    DiversionDiffSeg(:,:,1,MC)     =  MeanDivDiffLog;

    OneDiversionSameSeg(:,:,1,MC)     =  OneDivLog;
    OneDiversionDiffSeg(:,:,1,MC)     =  OneDivDiffLog;        
        
    FunctionValue(:,:,1,MC)         = FctValLog;
    ExitFlag(:,:,1,MC)              = nan;
    AIC(:,:,1,MC)                   = AICLogit;
    BIC(:,:,1,MC)                   = BICLogit;
    
    %% Nested Logit
    XNestedLog  = [Xexo lnsjg];
    BNLog       = ((XNestedLog'*Z*W*Z'*XNestedLog)\eye(size(XNestedLog,2)))*(XNestedLog'*Z*W*Z'*y);
    % Standard errors
    est         = y-XNestedLog*BNLog;
    dgf         = (size(XNestedLog,1)-size(Xexo,2));
    ser         = (est'*est)./dgf;
    sst         = inv(XNestedLog'*XNestedLog);
    sterrNL     = sqrt(ser*diag(sst));
    deltNL      = y - BNLog(nlin+1)*lnsjg;
    % Nested Logit elasticity wrt X1
    alpha       = BNLog(nlin);
    theta       = [0;0;BNLog(nlin+1)];
    [ElaNL, DiversionNL]     = ElastNestedLogit(alpha,theta,deltNL,multiX1,Data);
    
    [MeanOwnElastNL,OneOwnElastNL,MeanCrossElastSameSegNL,OneCrossElastSameSegNL,MeanCrossElastDiffSegNL,OneCrossElastDiffSegNL,...
     MeanDivNL,MeanDivDiffNL,OneDivNL,OneDivDiffNL] = sumElast(ElaNL,DiversionNL,Data);    
    % Function value
    FctValNL   = (est' * Z) * W * (Z' * est);
    AICNLog    = nobs * FctValNL - 2 * (size(Z,2)-size(BNLog,1));
    BICNLog    = nobs * FctValNL - (size(Z,2)-size(BNLog,1)) * log(nobs);
    
    % Store the results
    Coefficients(:,:,2,MC)          = [BNLog;nan;nan];
    StandardErrors(:,:,2,MC)        = [sterrNL;nan;nan];
    ComputationTime(:,:,2,MC)       = nan;
    OwnElast(:,:,2,MC)              = MeanOwnElastNL;
    CrossElastSameSeg(:,:,2,MC)     = MeanCrossElastSameSegNL;
    CrossElastDiffSeg(:,:,2,MC)     = MeanCrossElastDiffSegNL;
    
    OneOwnElast(:,:,2,MC)           = OneOwnElastNL;
    OneCrossElastSameSeg(:,:,2,MC)  = OneCrossElastSameSegNL;
    OneCrossElastDiffSeg(:,:,2,MC)  = OneCrossElastDiffSegNL;
          
    DiversionSameSeg(:,:,2,MC)     =  MeanDivNL;
    DiversionDiffSeg(:,:,2,MC)     =  MeanDivDiffNL;

    OneDiversionSameSeg(:,:,2,MC)     =  OneDivNL;
    OneDiversionDiffSeg(:,:,2,MC)     =  OneDivDiffNL;
        
    FunctionValue(:,:,2,MC)         = FctValNL;
    ExitFlag(:,:,2,MC)              = nan;
    AIC(:,:,2,MC)                   = AICNLog;
    BIC(:,:,2,MC)                   = BICNLog;
    
    %% Random Coefficient Logit I with random coefficient on continuous characteristic
    % Some data that speeds up computations
    xzwz            = Xexo'*Z*W*Z';
    Data.xzwz       = xzwz;
    xzwzx           = xzwz * Xexo;
    invxzwzx        = xzwzx\eye(size(xzwzx,2));
    Data.invxzwzx   = inv(xzwzx);
    Data.Z          = Z;
    Data.W          = W;
    Data.Xexo       = Xexo;
    % OPTIMIZATION
    % Options:see Knitro Option File KnitroBLP.opt
    % Pass Data into gmm using anonymous functions
    theta0RCL  = 1;     % starting value
    delta0     = deltL;
    save mvalold delta0
    
    angmm    = @(theta0RCL)gmmLogit(theta0RCL,Data);
    
    options  = optimset( 'Display','iter',...
        'GradObj','on','TolCon',1E-6,...
        'TolFun',1E-6,'TolX',1E-6,...
        'Hessian', 'off','DerivativeCheck','off','FinDiffType','central');
    
    t1      = cputime;
    [thetaRCL, FctValRCL, exitflagRCL, ~, ~] = ...
        ktrlink(angmm,theta0RCL, [], [], [], [], [], [], [], options, 'knitroBLP.opt');
    
    % Ricalculate beta according to the thetaRCL that minimizes the function
    % given the use of multistart by KNITRO
    
    deltRCL     = deltaLogit(thetaRCL,Data);
    beta        = invxzwzx*(xzwz*deltRCL);
    th12RCL     = [beta ; thetaRCL];
    sterrRCL    = seRCLogit(th12RCL,Data);
    cputimegmm  = cputime-t1;
    
    % Random Coefficient Logit elasticity wrt X1
    alpha        = beta(nlin);
    [ElaRCLogit,DiversionRCLogit] = ElastLogit(alpha,thetaRCL,deltRCL,multiX1,Data);
    
    [MeanOwnElastRCL,OneOwnElastRCL,MeanCrossElastSameSegRCL,OneCrossElastSameSegRCL,MeanCrossElastDiffSegRCL,OneCrossElastDiffSegRCL,...
     MeanDivRCL,MeanDivDiffRCL,OneDivRCL,OneDivDiffRCL] = sumElast(ElaRCLogit,DiversionRCLogit,Data);
    
    AICRCLog     = nobs * FctValRCL - 2 * (size(Z,2)-size(th12RCL,1));
    BICRCLog     = nobs * FctValRCL - (size(Z,2)-size(th12RCL,1)) * log(nobs);
    
    
    % Store the results
    Coefficients(:,:,3,MC)          = [th12RCL(1:nlin);nan;nan;abs(thetaRCL)];
    StandardErrors(:,:,3,MC)        = [sterrRCL(1:nlin);nan;nan;sterrRCL(nlin+1)];
    ComputationTime(:,:,3,MC)       = cputimegmm;
    OwnElast(:,:,3,MC)              = MeanOwnElastRCL;
    CrossElastSameSeg(:,:,3,MC)     = MeanCrossElastSameSegRCL;
    CrossElastDiffSeg(:,:,3,MC)     = MeanCrossElastDiffSegRCL;
       
    OneOwnElast(:,:,3,MC)           = OneOwnElastRCL;
    OneCrossElastSameSeg(:,:,3,MC)  = OneCrossElastSameSegRCL;
    OneCrossElastDiffSeg(:,:,3,MC)  = OneCrossElastDiffSegRCL;
             
    DiversionSameSeg(:,:,3,MC)     =  MeanDivRCL;
    DiversionDiffSeg(:,:,3,MC)     =  MeanDivDiffRCL;

    OneDiversionSameSeg(:,:,3,MC)     =  OneDivRCL;
    OneDiversionDiffSeg(:,:,3,MC)     =  OneDivDiffRCL;
        
    FunctionValue(:,:,3,MC)         = FctValRCL;
    ExitFlag(:,:,3,MC)              = exitflagRCL;
    AIC(:,:,3,MC)                   = AICRCLog;
    BIC(:,:,3,MC)                   = BICRCLog;
    
    %% Random Coefficient Logit II with random coefficient on continuous characteristic + discrete characteristic
    t1      = cputime;
    % OPTIMIZATION
    % Options:see Knitro Option File KnitroBLP.opt
    % Pass Data into gmm using anonymous functions
    theta20RCL     = [1;1];
    delta0         = deltL;
    save mvalold delta0
    x_L     = [-10,-10];         % lower bound to avoid numerical problems when sigma gets negative: negative weights make the shares negative is their single value is too small
    x_U     = [10,10];       % upper bound is really not binding

    angmm    = @(theta20RCL)gmmLogit(theta20RCL,Data);
    
    options  = optimset( 'Display','iter',...
        'GradObj','on','TolCon',1E-6,...
        'TolFun',1E-6,'TolX',1E-6,...
        'Hessian', 'off','DerivativeCheck','off','FinDiffType','central');
    
%     [thetaRCL2, fval, exitflag, output, lambda] = ...
%     fminunc(angmm,theta20RCL,options);
    
    [thetaRCL2, FctValRCL2, exitflagRCL2, ~, ~] = ...
        ktrlink(angmm,theta20RCL, [], [], [], [], x_L, x_U, [], options, 'knitroBLP.opt');
    
    % Ricalculate beta according to the thetaRCL that minimizes the function
    % given the use of multistart by KNITRO
    
    deltRCL2      = deltaLogit(thetaRCL2,Data);
    beta2         = invxzwzx*(xzwz*deltRCL2);
    th12RCL2      = [beta2 ; thetaRCL2];
    sterrRCL2     = seRCLogit(th12RCL2,Data);
    cputimegmm    = cputime-t1;
    
    % Random Coefficient Logit elasticity wrt X1
    alpha2        = beta2(nlin);
    [ElaRCLogit2,DiversionRCLogit2] = ElastLogit(alpha2,thetaRCL2,deltRCL2,multiX1,Data);
    
    [MeanOwnElastRCL2,OneOwnElastRCL2,MeanCrossElastSameSegRCL2,OneCrossElastSameSegRCL2,MeanCrossElastDiffSegRCL2,OneCrossElastDiffSegRCL2,...
     MeanDivRCL2,MeanDivDiffRCL2,OneDivRCL2,OneDivDiffRCL2] = sumElast(ElaRCLogit2,DiversionRCLogit2,Data);
 
    AICRCLog2     = nobs * FctValRCL2 - 2 * (size(Z,2)-size(th12RCL2,1));
    BICRCLog2     = nobs * FctValRCL2 - (size(Z,2)-size(th12RCL2,1)) * log(nobs);
    
     % Store the results
    Coefficients(:,:,4,MC)          = [th12RCL2(1:nlin);nan;abs(thetaRCL2)];
    StandardErrors(:,:,4,MC)        = [sterrRCL2(1:nlin);nan;sterrRCL2(nlin+1);sterrRCL2(nlin+2)];
    ComputationTime(:,:,4,MC)       = cputimegmm;
    OwnElast(:,:,4,MC)              = MeanOwnElastRCL2;
    CrossElastSameSeg(:,:,4,MC)     = MeanCrossElastSameSegRCL2;
    CrossElastDiffSeg(:,:,4,MC)     = MeanCrossElastDiffSegRCL2;
       
    OneOwnElast(:,:,4,MC)           = OneOwnElastRCL2;
    OneCrossElastSameSeg(:,:,4,MC)  = OneCrossElastSameSegRCL2;
    OneCrossElastDiffSeg(:,:,4,MC)  = OneCrossElastDiffSegRCL2;
             
    DiversionSameSeg(:,:,4,MC)     =  MeanDivRCL2;
    DiversionDiffSeg(:,:,4,MC)     =  MeanDivDiffRCL2;

    OneDiversionSameSeg(:,:,4,MC)     =  OneDivRCL2;
    OneDiversionDiffSeg(:,:,4,MC)     =  OneDivDiffRCL2;
        
    FunctionValue(:,:,4,MC)         = FctValRCL2;
    ExitFlag(:,:,4,MC)              = exitflagRCL2;
    AIC(:,:,4,MC)                   = AICRCLog2;
    BIC(:,:,4,MC)                   = BICRCLog2;
    
    
   %% Random Coefficient Nested Logit: this is the correctly specified DGP
    
    t1      = cputime;
    
    % OPTIMIZATION
    % Options:see Knitro Option File KnitroBLP.opt
    % Pass Data into gmm using anonymous functions
    theta20RCNL  = [0.5;0.5];     % starting values - just pick a low value of pho to avoid overflow problems from the start
    delta0NL     = deltaTrue;
    save mvaloldNL delta0NL
    
    angmmNL = @(theta20RCNL)gmmNL(theta20RCNL,Data);
    
    options = optimset( 'Display','iter',...
        'GradObj','on','TolCon',1E-6,...
        'TolFun',1E-6,'TolX',1E-6,...
        'Hessian', 'off','DerivativeCheck','off','FinDiffType','central');
    
    %     [thetaRCNL, FctValRCNL, exitflagRCNL, output, gradient] = ...
    %         fminunc(angmmNL,theta20RCNL,options);
    
    x_L     = [0;0];           % lower bound (only for the nested logit, to avoid numerical problems when pho gets negative)
    x_U     = [10;0.9];       % upper bound (only for the nested logit, to avoid numerical problems when pho gets close to 1)
    [thetaRCNL, FctValRCNL, exitflagRCNL,~,~] = ...
        ktrlink(angmmNL,theta20RCNL, [], [], [], [], x_L, x_U, [], options, 'knitroBLP.opt');
    
    deltRCNL      = deltaNL(thetaRCNL,Data);
    betaNL        = invxzwzx*(xzwz*deltRCNL);
    th12RCNL      = [betaNL ; abs(thetaRCNL)];
    sterrRCNL     = seNL(th12RCNL,Data);
    cputimegmm    = cputime-t1;
    
    % Random Coefficient Nested Logit wrt X1
    alpha         = betaNL(nlin);
    
    [ElaRCNL,DiversionRCNL] = ElastNestedLogit(alpha,thetaRCNL,deltRCNL,multiX1,Data);
    
    [MeanOwnElastRCNL,OneOwnElastRCNL,MeanCrossElastSameSegRCNL,OneCrossElastSameSegRCNL,MeanCrossElastDiffSegRCNL,OneCrossElastDiffSegRCNL,...
     MeanDivRCNL,MeanDivDiffRCNL,OneDivRCNL,OneDivDiffRCNL] = sumElast(ElaRCNL,DiversionRCNL,Data);
 
    AICRCNLog    = nobs * FctValRCNL - 2 * (size(Z,2)-size(th12RCNL,1));
    BICRCNLog    = nobs * FctValRCNL - (size(Z,2)-size(th12RCNL,1)) * log(nobs);
    
    % Store the results
    Coefficients(:,:,5,MC)          = [betaNL;thetaRCNL(2);nan;thetaRCNL(1)];
    StandardErrors(:,:,5,MC)        = [sterrRCNL(1:3);sterrRCNL(5);nan;sterrRCNL(4)];
    ComputationTime(:,:,5,MC)       = cputimegmm;
    OwnElast(:,:,5,MC)              = MeanOwnElastRCNL;
    CrossElastSameSeg(:,:,5,MC)     = MeanCrossElastSameSegRCNL;
    CrossElastDiffSeg(:,:,5,MC)     = MeanCrossElastDiffSegRCNL;
               
    OneOwnElast(:,:,5,MC)           = OneOwnElastRCNL;
    OneCrossElastSameSeg(:,:,5,MC)  = OneCrossElastSameSegRCNL;
    OneCrossElastDiffSeg(:,:,5,MC)  = OneCrossElastDiffSegRCNL;
    
    DiversionSameSeg(:,:,5,MC)     =  MeanDivRCNL;
    DiversionDiffSeg(:,:,5,MC)     =  MeanDivDiffRCNL;

    OneDiversionSameSeg(:,:,5,MC)     =  OneDivRCNL;
    OneDiversionDiffSeg(:,:,5,MC)     =  OneDivDiffRCNL;
        
    FunctionValue(:,:,5,MC)         = FctValRCNL;
    ExitFlag(:,:,5,MC)              = exitflagRCNL;
    AIC(:,:,5,MC)                   = AICRCNLog;
    BIC(:,:,5,MC)                   = BICRCNLog;

    
    
   savefile = '7resultsRCNL_HighCorr_Pho03Sigma1_Dstar12RandomCoeff.mat';
   save    (savefile, 'Coefftrue', 'Coefficients', 'StandardErrors', 'ComputationTime',...
        'OwnElast', 'CrossElastSameSeg', 'CrossElastDiffSeg',...
        'OneOwnElast','OneCrossElastSameSeg','OneCrossElastDiffSeg',...
        'DiversionSameSeg', 'DiversionDiffSeg', 'OneDiversionSameSeg','OneDiversionDiffSeg',...
        'BiasOwnElast', 'BiasCrossElastSameSeg', 'BiasCrossElastDiffSeg',...
        'FunctionValue', 'ExitFlag', 'AIC', 'BIC', ...
        'SaveSegment', 'SaveX1', 'OutsideShare', 'corrDstarX1')
end